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Abstract 

In this paper we describe splitting methods for solving Levitron, which 
is motivated to simulate magnetostatic traps of neutral atoms or ion traps. 
The idea is to levitate a magnetic spinning top in the air repelled by a 
base magnet. 

The main problem is the stability of the reduced Hamiltonian, while 
it is not defined at the relative equilibrium. Here it is important to derive 
stable numerical schemes with high accuracy. For the numerical studies, 
we propose novel splitting schemes and analyze their behavior. We deal 
with a Verlet integrator and improve its accuracy with iterative and ex- 
trapolation ideas. Such a Hamiltonian splitting method, can be seen as 
geometric integrator and saves computational time while decoupling the 
full equation system. 

Experiments based on the Levitron model are discussed. 

Keywords splitting method, Verlet integrator, iterative and extrapolation 
methods, Levitron problem. 

AMS subject classifications. 65M12, 65L06, 65P10. 

1 Introduction 

We are motivated to simulate a Levitron, which is a magnetic spinning top and 
can levitate in a magnetic field. The main problem of such a nonlinear problem 
is to achieve a stability for the calculation of the critical splint rate. While the 
stability of Levitrons are discussed in the work of [3] and their dynamics in [?], 
we concentrate on improving the standard time-integrator schemes for the re- 
duced Hamiltonian systems. It is important to derive stable numerical schemes 
with high accuracy to compute the non-dissipative equation of motions. For the 
numerical studies, we propose novel splitting schemes and analyze their behav- 
ior. We deal with a standard Verlet integrator and improve its accuracy with 
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iterative and extrapolation ideas. Such a Hamiltonian splitting method, can be 
seen as geometric integrator and saves computational time while decoupling the 
full equation system, see the splitting ideas in the overview article pQ. 

In the following we describe the reduced model of Gans [5] and an extension 
based on a novel idea of magentic field of Dullin [3 for a disk. 

1.1 Hamiltonian of Gans 

In the paper, we deal with the following problem (reduced Hamiltonian): 

1 / 2 2 2 pi (p 5 - P - 6cosg 4 ) 2 p 6 



I \ a asm oa c 



M 



a sin <74 



sin?4 cos q 5 — sin q 5 — + cosg 4 — 
dqi dq 2 J dq 3 



<?3 (1) 



The evolution of the dynamical variable u(q, p) (including q and p them- 
selves) is given by the Poisson bracket, 

n , . f du dH du dH\ 

^ (q ' P) -(aq-ap-ap-aq) =(A + Wq ' P) - (2) 
For the non-separable Hamiltonian of ([I]), we have: 

<9i/ . _ / p4 ( P 5-p 6 cos g 4 ) 2 P6 (cos 2 g 4 + (a/c) sin|) -p 5 cos q4 \ ,„n 

IP.qj- 1P1 ; P2,P3, „ , asin 2 94 > asin 2 g 4 



<9p 



The same is given for: 



dH 

5q" (P ' q) 

/ g2y o 2 y \ „ / . a 2 * a 2 * 

(M sin g 4 cos g 5 —-5- + cos g 4 q — ^— , M sm q A cos g 5 — -5- + cos g 4 - — — 

V 9gf oqiOq 3 J V 9< ?2 0(720q 3 

/ / a 2 * a 2 * \ a 2 *\ 

M ( sm q A [ sm g 5 a ^ a ^ + cos g 5 a _ a ^ ) + cos g 4 -^-g- ) - 1 



dq 2 dq 3 dqidq s 
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* . p 6 (p 5 -p 6 cosq 4 ) cosg 4 (p 5 -p 6 cosq 4 ) 2 

M cos q4 sm q 5 h cos q$ - — — sin q± - ' 



92 dqi J dq 3 J a sin g 4 a sin 3 g 4 

d$> d^\\ 

M sin g 4 cos q 5 — sing 5 — ,0) (4) 



dq 2 dqi 
A and B are Lie operators, or vector fields 

3p 9q ^ <9q <9p ^ 

The transfer to the operators are given in the following description. 
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The exponential operators e hA and e hB are then just shift operators, with 
Ti(}i) is a symmetric second order splitting method: 

T 2< w(h)(At) = e (A*/2)B e A t A e (A t /2)B 

and corresponds to the velocity form of the Verlet algorithm ( VV) . 
Further the splitting scheme: 

T 2 ,P V (h)(At) = e ( At /2)A e AtB e (A t /2)A (?) 

and corresponds to the position- form of the Verlet algorithm (PV) . 

See also the derivation of the Verlet algorithm in Appendix [7J 

Tiyv (h)(At) — Sab(Ii), the symplectic Verlet or leap-frog algorithm is given 

as: 

We start with (qcbPo)* = (q(* n ), pO"))': 

1 f)ff r) 

(q ljPl ) t =e' l / 2S (q ,p )* = (I - °— ( Pi , q^Xqo, Po)*, (8) 

(q 2 ,v 2 ) t =e^(q 1 ,v 1 ) i = (/ + |_(p, ; , q 4 )^-)(q 1; Pl ) 4 , (9) 



(q3,v 3 ) t =e^ (q2iP2)t = {I __ h J2^-(p i ,q i )^ i )(q 2 ,p 2 ) t . (10) 

And the substitution is given the algorithm for one time-step n — > n + 1 and 
we obtain: 

(q(t"+ 1 ),v(r+ 1 )) t = (q 3 ,V3) t . 

1.2 Higher order Expansion of Ver let-algorithm 

In the following we extend the Verlet algorithm with respect to higher order 
terms. 

Such terms are important in the application with iterative schemes to achieve 
higher order schemes. 

We start with (q ,Po)' = (q(* n )> p(*"))*: 



(qj.Pi)* = e' l / 2S (q ,p )* 
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(q 2 ,v 2 ) 4 = e^qx.vx)* 

= ( / + E7 ! (^Ea p (p-^)a q :) J )( ( io,Po) t , (12) 



(q 3 ,v 3 ) 4 = e h ^ B ( q2 ,p 2 r 

- (i + tk~l^^(p^)^y)( q(hPo y. (i3) 

And the substitution is given the algorithm for one time-step n — > n + 1 and 
we obtain: 

(q(i"+ 1 ),v(i" + 1 ))* = (q 3 ,V3) t . 



2 Iterative Schemes for coupled problems 

Based on the nonlinear equations, we have to deal with linearization or nonlinear 
averaging techniques. In the following, we discuss the fixed point iteration and 
Newton's method. 

We solve the nonlinear problem: 

F{x) = 0, (14) 

where F : W 1 ->■ W 1 . 



2.1 Fixed-point iteration 

The nonlinear equations can be formulated as fixed-point problems: 

x = K(x), (15) 

where K is the fixed-point map and is nonlinear, e.g. K(x) = x — F{x). 
A solution of ( 16 ) is called fix-point of the map K. 
The fix-point iteration is given as: 

= K(xi), (16) 

and is called nonlinear Richardson iteration, Picard iteration, or the method of 
successive substitution. 

Definition 2.1 Let il < K" and let G : Vt — > K m . G is Lipschitz continuous 
on fi with Lipschitz constant 7 if 

||G(a;)-G(i/)||<7||x-i/||, (17) 

for all x,y e fi. 
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For the convergence we have to assume that K be a contraction map on Q 
with Lipschitz constant 7 < 1. 

Algorithm 2.1 We apply the fix-point iterative scheme to decouple the non- 
separable Hamiltonian problem t3V and Op. 



qi=g-(Pi-i,qi-i),i6[i",(" +1 ] (18) 

a it 

Pi = -g-(Pi-l,q i -l),tG[t n ,t B+1 ] (19) 
p(t )=p ,q(t°)=q , (20) 

the starting solutions for the i-th iterative steps are given as : 
Pi—i(t), qj_i(t) are the solutions of the i — 1 th iterative step 
and we have the initial condition for the fix-point iteration: 

(po(t),qo(i)) t = (p(i n ),q(* n ))' 

We assume that we have convergent results after i = 1 . . . , m iterative steps 
or with the stopping criterion: 

max(||p i+ i - pi||, ||q i+1 - q^l) < err, 

while || • || is the Euclidean norm (or a simple vector-norm, e.g. L 2 ). 

Iterative Verlet applied to the Hamiltonian (J3J) and Q : 
We start with (q ,Po)' = (q(* n )> p(* n ))*: 
The iterative scheme is given as: 



q*(t) = q(0 + ^(p(t")-^|-(p i -i(t),q i -i(t)),q(* n )) 1 (21) 
h dH 

Pi(t) = p(t n )--g-(p i _i(t),q i _i(t)) 

hdH ( ( . „, hdH , . , . ,.\ 
-2^(K)-2aq- (Pi - l(t) ' qi - l(t)) j' 

q(n + ^ ( p(n ~ 2^ (p *- l(i); *-i(*)).q(* n ))^, (22) 

forte [t n ,t n+1 ],/i = i n+1 -t n ,ra = 0,l,...,iV, (23) 
1 = 1,2,3...,/, (24) 

where / = 3 or 4. 



For the fix-point iteration, we have the problem of the initialization, means 
the start of the iterative scheme. We can improve the starting solution with a 
preprocessing method, which derives a first improved initial solution. 

Improve Initialization Process 
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Algorithm 2.2 To improve the initial solution we can start with: 

1. ) We initialize with a result of the explicit Euler-method : 

(QOjPo)* = (fi{t n+1 ) EulerlstiP{t n+1 ) EulerlstY ■' 

2. ) We initialize with a result of the explicit RK-method : 
(qo,Po)* = (q(t n+1 )w 4t h,p(i" +1 )^4 t / l ) t .- 

2.2 Newton's method 



We solve the nonlinear operator equation ( 14 ) . 

While F:DcI->y with the Banach spaces X, Y is given with the norms 
|| • \ \x and || • \\y- Let F be at least once continuous differentiable, further we 
assume xq is a starting solution of the unknown solution x*. 

Then the successive linearization lead to the general Newton's method: 

F'{ Xi )Axi = -F{ Xi ), (25) 

where Axi = Xi+± — Xi and i = 0, 1, 2, 

The method derive the solution of a nonlinear problem by solving the fol- 
lowing algorithm. 

Algorithm 2.3 By considering the sequential splitting method we obtain the 
following algorithm. We apply the equations 

q-^=-(p,q)=0,t€ [t n ,t n+1 ] (26) 

p + !_(p,q)=0,te [t n ,t n+1 ] (27) 
P„ = p(* n ),qn = q(n (28) 
where p = (pi, . . . an d q = . . . into the Newtons- formula we have: 

F(p,q)=x+| x (x) (29) 

and we can compute 

x (fc+i) = x (fe) _ J D(F(x( fe )))^ 1 F(x< fe )), (30) 

where D(F(x.)) is the Jacobian matrix and k = 0, 1, . . .. 

We stop the iterations when we o Main : |x( fe+1 ) - xW| < err, where err is 
an error bound, e.g. err = 10~ 4 . 



3 SPLITTING METHODS 



The solution vector F is given as: 



F(x) 



/ ^,i(x) \ 
^, 2 (x) 

^,e(x) 
F p , 2 ( x ) 

V ^P,6(x) J 



where x = (q u . . . , g 6 ,Pi, • • • ,PeY and 



1 dp 



The Jacobian matrix for the equation system is given as 



dF\ 




dF\ 


dx 1 


X 2 


dx 12 


dF 2 


dF 2 


0F 2 


dx± 


X 2 


dx 12 


dF 12 


dF 12 


dF 12 


dxi 


x 2 


dx 12 



L>F(x) 



V 

where x= (x 1 , . . . , x Y2 f = (q x , . . . , q 6 , p 1: . . . ,p 6 ) f . 

3 Splitting Methods 

In the following, we discuss the different splitting schemes. 
The simplest such symmetric product is 

T 2 (h) = S AB (h) or T 2 (h) = S BA (h). 
If one naively assumes that 

T 2 (h) = c At ^ A+B '> +Ch 3 + Dh 4 + --- , 
then a Richardson extrapolation would only give 

^ [er 2 k (h/k) T 2 (h)] = e A *^+ B ) + o(h% 

a third-order algorithm. 
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Thus for a given set of n whole numbers {fci} one can have a 2nth-order 
approximation 

e At(A+B) = J- CO** + 0(h 2n+1 ). (37) 

For orders four to ten, one has explicitly: 

%(At) = - 1 -T 2 (At) + ^T 2 2 (^y (38) 

1 o-/a^ , 16 -r2 / A A 729 3 ^ At ^ , 1024^ {At\ 



7S(A t ) = --75(A t ) + -7? j - ^7? ^-j- j + ^7? j , (40) 

16384 4 /AA 390625 5 f At\ 
_ ^835" T2 l^Tj+^Ty 72 l^Tj' ( lj 

4 Iterative MPE method 

Based on the nonlinear problem, we extend the MPE method to an iterative 
scheme. 

Algorithm 4.1 Iterative Verlet applied to the Hamiltonian and : 
We start with (q ,Po)* = (q(* n ), p(* n ))*-" 

The iterative scheme for computing TiiPiiQi-, h) — (qi(/i), p$(fa) ) is given as: 



nPi-V„h) = e<" 2B e" e '" 2 ' J (q„,p )' 

1 9 



3=1 

JY 



j\ v dp ^ l - 1 ^ l - 1 'dq(t- 
1 / 1 . dH , ,9 



( J + Ej!(-2^^- 1 '*- 1 )9 P T^) ? )^ n )' p ^ 
VKe /lave i/ie higher order schemes given as: 
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%{pi,quh) = -~T2(pi,qi,h) + ^T 2 2 (pi^i, ~\ (43) 

m i \ ,n 16^/ 'A 81 / h\ 

T 6 [pi,qi,h) = ^T 2 (Pi,qi,h) - —T 2 [P"1i> 2 ) + 40 2 { Pi,qi, 3J ' 



. . 1 . . 16 2 / ^ 

Ts(Pi,qi,h) = -—T 2 {p i ,q l ,h) + —T 2 \p l -,q l , - 

729^o / /A 1024 . / h\ 
~ WQ Ti f Rlftl - +— tMp^-) (45) 



/ 1 o-/ ,n 64 „ / /A 6561 / ft 

71o(Pi,ft,/»)=8g S Qr a (Pi ) « i ,/»)-^ ^^-2j + 4480 r2 (, M '3 

16384 . / ft\ 390625^./ h\ 
"W 75 ( v Pl '*'4j + T2576 rT2 ^*'5j- (46) 

where forte [t n ,t n+1 ], h = t n+1 - t n ,n = 0,1, ... ,N, 
1 = 1,2,3...,/, 

and we /lave a stopping criterion or a fixed number of iterative steps, e.g. 
I = 3 or 4. 



5 Numerical Examples 

The Levitron is described on the base of rigid body theory. With the convention 
of Goldstein [5] for the Euler angles the angular velocity ujq is along the z-axis 
of the system, ojg along the line of nodes and ui^ along the z'-axis. Transforming 
them into body coordinates one gets 

( 4> sin 9 sin ip + 9 cos ifA 
> sin 9 cos ip + 9 sin tp | (47) 
) cos 9 + tp 

Finally the kinetic energy can be written as 
1 



T = 
2 



m[x 2 + y 2 + z 2 ) + A{9 2 + <p 2 sin 2 9) + C(ip + <j> cos 9) 2 (48) 



The potential energy U is given by the sum of the gravitational energy and 
the interaction potential of the Levitron in the magnetic field of the base plate: 

$ $ $ 

U = mgz — //(sin tp sin 6 |-cosi/>sin0 hcos# — ) (49) 

x y z 

with mu as the magnetic moment of the top and $ the magneto-static poten- 
tial. Following Gans [6] we uses the potential of a ring dipole as approximation 
for a magnetized plane with a centered unmagnetized hole. Furthermore we 
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introduced a nondimensionalization for the variables and the magneto-static 
potential: 



* = 



-(X 2 



3 (2Z 2 - 3)Z 



(50) 



(1 + £2)3/2 V M(l + Z 2 ) 7 /2 

Lengths were scaled by the radius R of the base plane, mass were measured 
in units of m and energy in units of mgh. Therefore the one time unit is y/R/g. 
So the dimensionless Hamiltonian with q = (X, Y, Z, 9, rp, 4>) is given by 



H 



Pi 



M 



vl 



sm<74 



9 

P3 



pi 

a 



(p 5 - P - 6cosg 4 ) 2 



cos <75^— sin 95 — 
oqi aq 2 



a sin <74 

cos g 4 



- qz 



(51) 
(52) 



with a and c as the nondimensionalized inertial parameter and M as the ratio 
of gravitational and magnetic energy. 

Solving the equations of motion ^ and Q numerically with the methods 
described above, one is able to plot the movement of the center of mass like it 
is shown in figure [T] In the plot the first 5 seconds of a stable trajectory were 
plotted. With a longer calculation we tested, that the top would levitate for 
more than one minute. 




Figure 1: stable trajectory of the center of mass in 3D and 2D 

The axis in the plot show the nondimensional variables X,Y and Z. The 
trajectory starts at the equilibrium point (qi, 92, 93) = (0,0,1.72). This tra- 
jectory was calculated with the fourth-order RungeKutta method with a small 
timestep of 10 -5 units of time, where one time unit is about 59,5ms, because 
of the nondimensionalization. For further considerations this trajectory will be 
used as a reference solution. 

The errors of the different time-steps with Runge-Kutta are given in Figure 

m 

The same equations were solved with the iterative Verlet algorithm de- 
scribed before. Due to the long computation time needed, we simulated only 
1000 timesteps and compare the trajectory with the reference solution from the 
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Figure 2: Errors of the numerical scheme: Runge-Kutta method (explicit 4th 
order) . 

Runge-Kutta algorithm. In figure [3] is shown how the trajectory of the same 
initial conditions looks like with the Verlet algorithm. 




-0.01 -0.008-0.006-0.004-0.002 0.002 0.004 0.006 0.008 



X 

Figure 3: trajectory calculated with Verlet algorithm 

This were done for one, two and four iterations per timestep, to see whether 
how many iterations are reasonable. The results are shown in Figure |4j 

In a first comparison, we deal with the second order Verlet algorithm and 
improve the scheme with iterative steps. 

In a first initialisation process, one can see that the errors are very similar 
to 1 or 2 iterative steps. 

The reduction of the error is possible with the improvement to higher order 
initialisation scheme, e.g. start with a first approximate solution with a RK 
scheme, or apply extrapolation schemes. 

The following tables [l] and [2] should give an impression of the timescales of 
the problem and the errors. 

Remark 5.1 Obviously the iterations does not improve the algorithm, when 
only using a lower order initialisation. By the way, it is sufficient to apply one 
iterative step in in comparison with the Runge-Kutta algorithm. 
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Figure 4: Errors of the numerical scheme: Iterative Verlet method. 



Runge-Kutta 


timestep 


1CT 5 


10- 3 


10" 1 


number of steps 


1000000000 


10000000 


100000 


computing time 


119 min 


23 sec 


2 sec 


stability 


ok 


ok 


ok 



Table 1: Stability and Computational Time with 4th order explicit Runge- 
Kutta. 



Also we have a benefit in reducing the computational time instead of applying 
only Runge-Kutta schemes. 

We tried to improve the solution with a extrapolation scheme in fourth 
order. We have a view at the errors this algorithm produces in comparison with 
the Runge-Kutta Solution with small time-steps (10 -5 time units per step). 
In Figure [5] we presented the results of the 4th MPE method with different 
time-steps and compared it with the Runge-Kutta solution. 





''.$1 


iUi ['■'.] :|? J-H 





Figure 5: Errors of the numerical scheme: 4th oder Extrapolation Scheme with 
Verlet method a Kernel (h = 10~ 5 ). 



Also we tested the 6th order MPE method with different time-steps and 
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Verlet 


timestep 


lCT 6 


10~ e 


10~ 6 


iterations per step 


1 


2 


4 


stability 


ok 


ok 


ok 


computing time 


67min 


120min 


219min 


mean error 


0.068 


0.068 


0.068 


maximal error 


0.0187 


0.0188 


0.0187 



Table 2: Stability and Computational Time with 2nd order Verlet Scheme. 



compared it with the Runge-Kutta solution, see Figure [6] 



Figure 6: Errors of the numerical scheme: 6th order Extrapolation Scheme with 
Verlet method a Kernel (h = 10~ 5 ). 



Like for Runge-Kutta we want to give an impression of the time scales for 
this extrapolation schemes, see Table [3] 





Extrapolation 4th order 


Extrapolation 6th order 


timestep 


10" 5 


10" e 


10" 5 


10" e 


number of steps 


100000000 


1000000000 


100000000 


1000000000 


computing time 


14min 


142min 


29min 


272min 


mean error 


0.007 


0.007 


0.0068 


0.0068 


maximal error 


0.0226 


0.0234 


0.0188 


0.0188 



Table 3: Errors and Computational Time with 4th order MPE scheme unsing 
Verlet Scheme as Kernel. 



6 Conclusions and Discussions 

In the paper, we have presented a model to simulate a Levitron. Based on 
the given Hamiltonian system, which is nonlinear, we present novel and simpler 
schemes based on splitting ideas to solve the equation systems. In future, we 
concentrate on the numerical analysis and embedding higher order splitting 
kernels to the extraopoation schemes. 
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7 Appendix 

For example, the evolution of any dynamical variable w(q, p) (including q and 
p themselves) is given by the Poisson bracket, 

dtU ^ p) = ( ai ' ap" - dp ' aq") = {A + B)u ^ p) - (53) 



For a separable Hamiltonian, 



ff(p,q) = ^+V(q). (54) 
A and £? are Lie operators, or vector fields 

A = v-A fl = a(q).A (55) 

where we have abbreviated §^-(p, q) = v = p/m and — ^(p,q) = a(q) = 
— W(q)/m. The exponential operators e 71 " 4 and e^ 5 are then just shift opera- 
tors. 

S{h) = c h/2B c hA c h/2B 

That is also given as a Verlet-algorithm in the following scheme. 
We start with (q ,v )* = (q(t rl ), v(t n ))*: 

(q 1 ,v 1 ) t =c' l / 2B (qo,vo)* = (/ + \h £ a(q)^-)(q , v )* (56) 

i 

= (qo,v + ^/ia(q )) t , (57) 



(q 2 ,v 2 ) t =c' l ' 4 (q 1 ,v 1 ) t = (/ + /i^v i ^-)(q 1 ,v 1 ) t (58) 

cq 4 

= (qi + ftvi.vi)*, (59) 



(q 3 ,v 3 ) t =c' l / 2B (q 2 ,v 2 )* = (/+^^a(q)A )(q2;V2) * (60 ) 

= (q 2 ,v 2 + ^Mqi)) t - (61) 

And the substitution is given the algorithm for one time-step n — > n + 1: 

h 2 h h h 

(q3, v 3 )* = (q + frv + y a(q ), v + -a(q ) + ^ a (qo + ^v + -a(q )))*, (62) 
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while (q(t"+ 1 ),v(^+ 1 ))* = (q3,V3) 4 . 
Iterative Verlet Algorithm 

In the abstract version of f|(Pi-i, qi-i?), -f|(p»-i,qi-i). 

Algorithm 7.1 We have the iterative Verlet Algorithm: 

1.) We start with the initialisation : (p (t n+1 ),q (t n + 1 )) t = {p(t n ), q(t"))' 
and i = 

2.,) TTie iterative step is given as: i = i + 1 and we /icwe: 

< +1 = q" + ^(p" - 2^(p?-i. Cli 1 ). qT-V). (63) 



P,(^ +1 )-p"-^(pr- + 1 1 I qr +1 ), (64) 
p^ n+1 ) - p" - ^(pJ^.q" + ^(p b - ^(pJ^.q^ 1 ).^ 1 )). 

w/iere /i = — t n is the local time-step. 
We compute the stopping criterion: 

max(||p(i™ +1 - p^ll, ||q(i™ +1 - q^ll) < err or we stop after i= I, while 
I is the maximal iterative step. 
3.) The result is given as: 
(p(t"+ 1 ),q(t"+ 1 )j t = ( Pi (i"+ 1 ),q i (t"+ 1 ))* 

and n = n+ 1 if n > N, while N is the maximal time-step, we stop 
else we go to step 1.) 
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